##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(RColorBrewer)
library(ggpubr)
library(scales)

##########################################################################################

option_list <- list(
    make_option(c("--njmu_file"), type = "character") ,
    make_option(c("--tcga_file"), type = "character") ,
    make_option(c("--utokyo_file"), type = "character") ,
    make_option(c("--tmucih_file"), type = "character") ,
    make_option(c("--tmucih_msi_file"), type = "character") ,
    make_option(c("--tcga_clinic_file"), type = "character") ,
    make_option(c("--oncoSG_msi_file"), type = "character") ,
    make_option(c("--tcga_msi_file"), type = "character") ,
    make_option(c("--oncoSG_clinic_file"), type = "character") ,
    make_option(c("--oncoSG_file"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    njmu_file <- "~/20220915_gastric_multiple/dna_combinePublic/baseTable/STAD_Info.addBurden.MSI_MSS.tsv"
    tcga_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/TCGA/TCGA_STAD.TMB.tsv"
    oncoSG_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/OncoSG/OncoSG_STAD.TMB.tsv"

    ## utokyo
    utokyo_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/utokyo/utokyo_STAD.TMB.tsv"

    ## TMUCIH
    tmucih_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/TMUCIH/TMUCIH_STAD.TMB.tsv"    

    ## MSI
    oncoSG_msi_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/OncoSG/OncoSG_STAD_msi.TMB.tsv"
    tcga_msi_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/TCGA/TCGA_STAD_msi.TMB.tsv"

    tcga_clinic_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/TCGA/clinical_PANCAN_patient_with_followup.tsv"
    oncoSG_clinic_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/OncoSG/OncoSG_STAD.followup.tsv"

    tmucih_msi_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/TMUCIH/TMUCIH_STAD_msi.TMB.tsv"    

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

njmu_file <- opt$njmu_file
tcga_file <- opt$tcga_file
oncoSG_file <- opt$oncoSG_file
tcga_msi_file <- opt$tcga_msi_file
tcga_clinic_file <- opt$tcga_clinic_file
oncoSG_clinic_file <- opt$oncoSG_clinic_file
oncoSG_msi_file <- opt$oncoSG_msi_file
out_path <- opt$out_path
utokyo_file <- opt$utokyo_file
tmucih_file <- opt$tmucih_file
tmucih_msi_file <- opt$tmucih_msi_file

dir.create(out_path , recursive = T)

###########################################################################################

if(1!=1){
    col <- c(
      brewer.pal(9,"YlGnBu")[6],
      rgb(234,106,79,alpha=255,maxColorValue=255),
      rgb(203,24,30,alpha=255,maxColorValue=255),
      rgb(255,0,0,alpha=255,maxColorValue=255)
      )
}
col <- c(
    brewer.pal(9,"YlGnBu")[6],
    rgb(red=247,green=184,blue=71,alpha=255,max=255) ,
    rgb(red=2,green=100,blue=190,alpha=255,max=255) ,
    rgb(255,0,0,alpha=255,maxColorValue=255)
    )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")
#col <- c("#E41A1C" , "#4DAF4A","#377EB8" ,"#984EA3")

###########################################################################################

dat_njmu <- data.frame(fread(njmu_file , quote=""))

dat_tcga <- data.frame(fread(tcga_file , quote=""))
dat_tcga_msi <- data.frame(fread(tcga_msi_file , quote=""))
dat_tcga_clnic <- data.frame(fread(tcga_clinic_file , quote=""))

dat_oncosg <- data.frame(fread(oncoSG_file , quote=""))
dat_oncosg_msi <- data.frame(fread(oncoSG_msi_file , quote=""))
dat_oncosg_clinic <- data.frame(fread(oncoSG_clinic_file , quote=""))

dat_utokyo <- data.frame(fread(utokyo_file , quote=""))

dat_tmucih <- data.frame(fread(tmucih_file , quote=""))
dat_tmucih_msi <- data.frame(fread(tmucih_msi_file , quote=""))


#dat_hk <- data.frame(fread(hk_file , quote=""))
#dat_hk_clinic <- data.frame(fread(hk_clinic_file , quote=""))

col_names <- c("Tumor" , "Age" , "Gender" , 
    "Tobacco" , "Alcohol" , "py" , "HP" ,
    "Stage" , "Class" , "From" , "MS_Type" , "BurdenExon")

###########################################################################################
## 去除同一个患者多个样本 
## 这5个样本保留
dat_njmu2 <- dat_njmu
dat_njmu <- subset( dat_njmu , Type!="IM + IGC + DGC")
dat_njmu <- dat_njmu[,c("Patient" , "Age" , "Gender" , 
    "Tobacco" , "Alcohol" , "py" , "HP" ,
    "Stage" , "Class" , "MS_Type" , "BurdenExon" )]
dat_njmu$From <- "NJMU"
dat_njmu <- subset( dat_njmu , Class %in% c("IGC" , "DGC") )
## 同一个人的合并
dat_njmu <- dat_njmu %>%
group_by( Patient , Age , Gender , Tobacco , Alcohol , py , HP , Stage , Class , From , MS_Type ) %>%
summarize( BurdenExon = median(BurdenExon) )
colnames(dat_njmu) <- col_names

###########################################################################################
## 同一患者均有的
dat_njmu2 <- subset( dat_njmu2 , Type=="IM + IGC + DGC")
dat_njmu2 <- subset( dat_njmu2 , Class %in% c("IGC" , "DGC") )
dat_njmu2$Class <- "IGC + DGC"
dat_njmu2 <- dat_njmu2[,c("Patient" , "Age" , "Gender" , 
    "Tobacco" , "Alcohol" , "py" , "HP" ,
    "Stage" , "Class" , "MS_Type" ,  "BurdenExon" )]
dat_njmu2$From <- "NJMU"
## 同一个人的合并
dat_njmu2 <- dat_njmu2 %>%
group_by( Patient , Age , Gender , Tobacco , Alcohol , py , HP , Stage , Class , From , MS_Type  ) %>%
summarize( BurdenExon = median(BurdenExon) )
colnames(dat_njmu2) <- col_names

dat_combine <- dat_njmu2
dat_combine$Gender[dat_combine$Gender %in% c("female" , "Female" , "FEMALE")] <- "Female"
dat_combine$Gender[dat_combine$Gender %in% c("male" , "Male" , "MALE")] <- "male"

dat_combine$Tobacco[is.na(dat_combine$Tobacco) | dat_combine$Tobacco == "" ] <- "unknown"
dat_combine$py[is.na(dat_combine$py) | dat_combine$py == ""] <- "unknown"
dat_combine$Alcohol[is.na(dat_combine$Alcohol) | dat_combine$Alcohol == "" ] <- "unknown"

dat_combine$Age[is.na(dat_combine$Age)] <- "unknown"
dat_combine$HP[is.na(dat_combine$HP)] <- "unknown"

dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("I" , "IA" , "IB" , "Stage IB" , "Stage IA") , "I" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("II" , "Stage II" , "Stage IIB" , "Stage IIA") , "II" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("III" , "IIIA" , "IIIB" , "Stage III" , "Stage IIIA" , "Stage IIIB", "Stage IIIC") , "III" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("Stage IV" , "IV") , "IV" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    !(dat_combine$Stage %in% c("I" , "II" , "III" , "IV")) , "unknown" , dat_combine$Stage
    )

dat_njmu2 <- dat_combine

###########################################################################################
## TCGA数据整理
dat_tcga$Class <- "DGC"
dat_tcga$Class[grep( "Intestinal" , dat_tcga$histological_type)] <- "IGC"
dat_tcga$BurdenExon <- dat_tcga$TMB..nonsynonymous.
dat_tcga$From <- "TCGA"

dat_tcga <- dat_tcga[,c("bcr_patient_barcode" , "gender" , "Class" , "From" , "BurdenExon")]

dat_tcga_clnic <- dat_tcga_clnic[,c(
    "bcr_patient_barcode" , "age_at_initial_pathologic_diagnosis" , "number_pack_years_smoked" ,
    "tobacco_smoking_history" , "alcohol_history_documented" , "pathologic_stage"
    )
]

dat_tcga_use <- merge( dat_tcga , dat_tcga_clnic )
dat_tcga_use$HP <- "unknown"
dat_tcga_use$MS_Type <- "MSS"
dat_tcga_use <- dat_tcga_use[,
    c("bcr_patient_barcode" , "age_at_initial_pathologic_diagnosis" , "gender" ,
    "tobacco_smoking_history" ,  "alcohol_history_documented" , "number_pack_years_smoked" , "HP" ,
    "pathologic_stage" , "Class" , "From" , "MS_Type" , "BurdenExon")
]

colnames(dat_tcga_use) <- col_names
dat_tcga_use_mss <- dat_tcga_use

##########################################
## TCGA数据MSI整理
dat_tcga <- dat_tcga_msi
dat_tcga$Class <- "DGC"
dat_tcga$Class[grep( "Intestinal" , dat_tcga$histological_type)] <- "IGC"
dat_tcga$BurdenExon <- dat_tcga$TMB..nonsynonymous.
dat_tcga$From <- "TCGA"

dat_tcga <- dat_tcga[,c("bcr_patient_barcode" , "gender" , "Class" , "From" , "BurdenExon")]

dat_tcga_clnic <- dat_tcga_clnic[,c(
    "bcr_patient_barcode" , "age_at_initial_pathologic_diagnosis" , "number_pack_years_smoked" ,
    "tobacco_smoking_history" , "alcohol_history_documented" , "pathologic_stage"
    )
]

dat_tcga_use <- merge( dat_tcga , dat_tcga_clnic )
dat_tcga_use$HP <- "unknown"
dat_tcga_use$MS_Type <- "MSI"
dat_tcga_use <- dat_tcga_use[,
    c("bcr_patient_barcode" , "age_at_initial_pathologic_diagnosis" , "gender" ,
    "tobacco_smoking_history" ,  "alcohol_history_documented" , "number_pack_years_smoked" , "HP" ,
    "pathologic_stage" , "Class" , "From" , "MS_Type" , "BurdenExon")
]

colnames(dat_tcga_use) <- col_names
dat_tcga_use_msi <- dat_tcga_use

##########################################

dat_tcga_use <- rbind(dat_tcga_use_msi , dat_tcga_use_mss)

###########################################################################################
## Oncosg
dat_oncosg$Class <- "DGC"
dat_oncosg$Class[grep( "Intestinal" , dat_oncosg$Laurens.classification)] <- "IGC"
dat_oncosg$BurdenExon <- dat_oncosg$TMB..nonsynonymous.
dat_oncosg$From <- "EastAsian"
dat_oncosg <- dat_oncosg[,c("Patient.ID" , "Sex" , "Class" , "From" , "BurdenExon")]

dat_oncosg_clinic <- dat_oncosg_clinic[,c(
    "Patient.ID" , "Age" , "Stage"
    )
]

dat_oncosg_use <- merge( dat_oncosg , dat_oncosg_clinic )
dat_oncosg_use$tobacco_smoking_history <- "unknown"
dat_oncosg_use$alcohol_history_documented <- "unknown"
dat_oncosg_use$number_pack_years_smoked <- "unknown"
dat_oncosg_use$HP <- "unknown"
dat_oncosg_use$MS_Type <- "MSS"

dat_oncosg_use <- dat_oncosg_use[,
    c("Patient.ID" , "Age" , "Sex" ,
    "tobacco_smoking_history" ,  "alcohol_history_documented" , "number_pack_years_smoked" , "HP" ,
    "Stage" , "Class" , "From" , "MS_Type" , "BurdenExon")
]

colnames(dat_oncosg_use) <- col_names
dat_oncosg_use_mss <- dat_oncosg_use

#######################################
## Oncosg的MSI
dat_oncosg <- dat_oncosg_msi
dat_oncosg$Class <- "DGC"
dat_oncosg$Class[grep( "Intestinal" , dat_oncosg$Laurens.classification)] <- "IGC"
dat_oncosg$BurdenExon <- dat_oncosg$TMB..nonsynonymous.
dat_oncosg$From <- "EastAsian"
dat_oncosg <- dat_oncosg[,c("Patient.ID" , "Sex" , "Class" , "From" , "BurdenExon")]

dat_oncosg_clinic <- dat_oncosg_clinic[,c(
    "Patient.ID" , "Age" , "Stage"
    )
]

dat_oncosg_use <- merge( dat_oncosg , dat_oncosg_clinic )
dat_oncosg_use$tobacco_smoking_history <- "unknown"
dat_oncosg_use$alcohol_history_documented <- "unknown"
dat_oncosg_use$number_pack_years_smoked <- "unknown"
dat_oncosg_use$HP <- "unknown"
dat_oncosg_use$MS_Type <- "MSI"

dat_oncosg_use <- dat_oncosg_use[,
    c("Patient.ID" , "Age" , "Sex" ,
    "tobacco_smoking_history" ,  "alcohol_history_documented" , "number_pack_years_smoked" , "HP" ,
    "Stage" , "Class" , "From" , "MS_Type" , "BurdenExon")
]

colnames(dat_oncosg_use) <- col_names
dat_oncosg_use_msi <- dat_oncosg_use

#######################################

dat_oncosg_use <- rbind( dat_oncosg_use_mss , dat_oncosg_use_msi )


###########################################################################################
## utokyo
dat_utokyo$Class <- "DGC"
dat_utokyo$BurdenExon <- dat_utokyo$TMB_NONSYNONYMOUS
dat_utokyo$From <- "Utokyo"
dat_utokyo <- dat_utokyo[,c("SAMPLE_ID" , "SEX" , "Class" , "From" , "BurdenExon" , "UICC_TUMOR_STAGE")]

dat_utokyo_use <- dat_utokyo
dat_utokyo_use$tobacco_smoking_history <- "unknown"
dat_utokyo_use$alcohol_history_documented <- "unknown"
dat_utokyo_use$number_pack_years_smoked <- "unknown"
dat_utokyo_use$HP <- "unknown"
dat_utokyo_use$MS_Type <- "MSS"
dat_utokyo_use$Age <- "unknown"

dat_utokyo_use <- dat_utokyo_use[,
    c("SAMPLE_ID" , "Age" , "SEX" ,
    "tobacco_smoking_history" ,  "alcohol_history_documented" , "number_pack_years_smoked" , "HP" ,
    "UICC_TUMOR_STAGE" , "Class" , "From" , "MS_Type" , "BurdenExon")
]

colnames(dat_utokyo_use) <- col_names
dat_utokyo_use[is.na(dat_utokyo_use)] <- "NA"


###########################################################################################
## TMUCIH
dat_tmucih$Class <- "DGC"
dat_tmucih$Class[grep( "I" , dat_tmucih$Subtype)] <- "IGC"
dat_tmucih$BurdenExon <- dat_tmucih$TMB..nonsynonymous.
dat_tmucih$From <- "TMUCIH"
dat_tmucih <- dat_tmucih[,c("Patient.ID" , "Diagnosis.Age" , "Stage" , "Sex" , "Class" , "From" , "BurdenExon")]

dat_tmucih_use <- dat_tmucih
dat_tmucih_use$tobacco_smoking_history <- "unknown"
dat_tmucih_use$alcohol_history_documented <- "unknown"
dat_tmucih_use$number_pack_years_smoked <- "unknown"
dat_tmucih_use$HP <- "unknown"
dat_tmucih_use$MS_Type <- "MSS"

dat_tmucih_use <- dat_tmucih_use[,
    c("Patient.ID" , "Diagnosis.Age" , "Sex" ,
    "tobacco_smoking_history" ,  "alcohol_history_documented" , "number_pack_years_smoked" , "HP" ,
    "Stage" , "Class" , "From" , "MS_Type" , "BurdenExon")
]

colnames(dat_tmucih_use) <- col_names
dat_tmucih_use_mss <- dat_tmucih_use

#######################################
## Oncosg的MSI
dat_tmucih <- dat_tmucih_msi
dat_tmucih$Class <- "DGC"
dat_tmucih$Class[grep( "I" , dat_tmucih$Subtype)] <- "IGC"
dat_tmucih$BurdenExon <- dat_tmucih$TMB..nonsynonymous.
dat_tmucih$From <- "TMUCIH"
dat_tmucih <- dat_tmucih[,c("Patient.ID" , "Diagnosis.Age" , "Stage" , "Sex" , "Class" , "From" , "BurdenExon")]

dat_tmucih_use <- dat_tmucih
dat_tmucih_use$tobacco_smoking_history <- "unknown"
dat_tmucih_use$alcohol_history_documented <- "unknown"
dat_tmucih_use$number_pack_years_smoked <- "unknown"
dat_tmucih_use$HP <- "unknown"
dat_tmucih_use$MS_Type <- "MSI"

dat_tmucih_use <- dat_tmucih_use[,
    c("Patient.ID" , "Diagnosis.Age" , "Sex" ,
    "tobacco_smoking_history" ,  "alcohol_history_documented" , "number_pack_years_smoked" , "HP" ,
    "Stage" , "Class" , "From" , "MS_Type" , "BurdenExon")
]

colnames(dat_tmucih_use) <- col_names
dat_tmucih_use_msi <- dat_tmucih_use

#######################################

dat_tmucih_use <- rbind( dat_tmucih_use_mss , dat_tmucih_use_msi )

###########################################################################################

#dat_combine <- rbind( data.frame(dat_njmu) , dat_oncosg_use , dat_tcga_use , dat_hk_use )
dat_combine <- rbind( data.frame(dat_njmu) , dat_oncosg_use , dat_tcga_use , dat_tmucih_use , dat_utokyo_use )

dat_combine$Gender[dat_combine$Gender %in% c("female" , "Female" , "FEMALE")] <- "female"
dat_combine$Gender[dat_combine$Gender %in% c("male" , "Male" , "MALE")] <- "male"
dat_combine$Gender[is.na(dat_combine$Gender) | dat_combine$Gender == "" | dat_combine$Gender == "NA"] <- "unknown"

dat_combine$Tobacco[is.na(dat_combine$Tobacco) | dat_combine$Tobacco == "" ] <- "unknown"
dat_combine$py[is.na(dat_combine$py) | dat_combine$py == ""] <- "unknown"
dat_combine$Alcohol[is.na(dat_combine$Alcohol) | dat_combine$Alcohol == "" ] <- "unknown"

dat_combine$Age[is.na(dat_combine$Age)] <- "unknown"
dat_combine$HP[is.na(dat_combine$HP)] <- "unknown"

dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("I" , "IA" , "IB" , "Stage IB" , "Stage IA" , "1") , "I" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("II" , "Stage II" , "Stage IIB" , "Stage IIA" , "2") , "II" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("III" , "IIIA" , "IIIB" , "Stage III" , "Stage IIIA" , "Stage IIIB", "Stage IIIC" , "3") , "III" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    dat_combine$Stage %in% c("Stage IV" , "IV" , "4") , "IV" , dat_combine$Stage
    )
dat_combine$Stage <- ifelse(
    !(dat_combine$Stage %in% c("I" , "II" , "III" , "IV")) , "unknown" , dat_combine$Stage
    )

###########################################################################################

dat_combine$From <- ifelse( dat_combine$From=="EastAsian" , "OncoSG" , dat_combine$From)

my_comparisons <- list( 
    c("OncoSG", "NJMU"), c("NJMU", "TCGA") , c("NJMU", "TMUCIH") , c("NJMU", "Utokyo") ,
    c("OncoSG", "TCGA") , c("OncoSG", "TMUCIH") , c("OncoSG", "Utokyo") ,
    c("TCGA", "TMUCIH") , c("TCGA", "Utokyo") ,
    c("TMUCIH", "Utokyo")
    )

dat_combine$From <- factor( dat_combine$From  , levels = c("NJMU" , "OncoSG" , "TCGA" , "TMUCIH" , "Utokyo" ) , order = T )
out_name <- paste0(out_path , "/MutationBurden.combine.tsv")
write.table( dat_combine , out_name , row.names = F , quote = F , sep = "\t" )

outcompare <- dat_combine %>%
group_by( Class , From , MS_Type ) %>%
summarize( BurdenExon = median(BurdenExon) )
out_name <- paste0(out_path , "/MutationBurden.combine.plot.tsv")
write.table( outcompare , out_name , row.names = F , quote = F , sep = "\t" )

###########################################################################################
## 包含了5个样本的
dat_combine2 <- rbind( dat_njmu2 , dat_combine)
out_name <- paste0(out_path , "/MutationInfo.combine.tsv")
write.table( dat_combine2 , out_name , row.names = F , quote = F , sep = "\t" )


###########################################################################################
## 比较不同来源的批次
y_lab <- "Mutation rate per MB"

for( msi in unique(unique(dat_combine$MS_Type)) ){

    dat_plot_tmp_use <- subset( dat_combine , MS_Type == msi )

    plot <- ggplot(dat_plot_tmp_use, aes(x = From , y = BurdenExon, color = From)) +
            geom_boxplot(size = 1.2 , outlier.alpha=0) + ## 去除散点，加粗线
            facet_grid(.~Class) + 
            geom_jitter(position = position_jitterdodge(0.2) , size = 1) + 
            scale_y_sqrt() +
            scale_color_npg() +
            xlab(NULL) +
            ylab(y_lab)+
            theme_bw() +
            #stat_compare_means(comparisons = my_comparisons , label = "p.format") +
            theme(
                legend.position = 'right',
                legend.title = element_blank() ,
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 8,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 10,color="black",face='bold',angle = 45, hjust = 0.5, vjust = 0.5) ,
                axis.line = element_line(size = 0.5)) 

    out_name <- paste0(out_path , "/MutationBurden.compare.From.",msi,".pdf")
    ggsave(file=out_name,plot=plot,width=8,height=6)
}

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("P == ",down," %*% 10","^",up)
    return(text)
}

###########################################################################################
## 分组做检验
for( msi in unique(unique(dat_combine$MS_Type)) ){

    dat_plot_tmp_use <- subset( dat_combine , MS_Type == msi )

    ###########################################################################################
    ## 样本数量
    sample_num <- dat_plot_tmp_use %>%
    group_by(From) %>%
    summarize( nums = length(unique(Tumor)) )

    sample_class_num <- dat_plot_tmp_use %>%
    group_by(Class) %>%
    summarize( nums_class = length(unique(Tumor)) )

    dat_plot_tmp_use <- merge( dat_plot_tmp_use , sample_num , by = "From")
    dat_plot_tmp_use <- merge( dat_plot_tmp_use , sample_class_num , by = "Class")


    dat_tmp <- c()
    for( from in unique(dat_plot_tmp_use$From) ){
        dat <- subset( dat_plot_tmp_use , From == from )
        if(from!="Utokyo"){
            a <- dat[dat$Class==unique(dat$Class)[1],"BurdenExon"]
            b <- dat[dat$Class==unique(dat$Class)[2],"BurdenExon"]
            
            p <- wilcox.test( a , b )$p.value

            if( p < 0.01 ){
                p_text <- trans(p)
            }else{
                p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
            }
        }else if(from=="Utokyo"){
            p_text <- ""
        }
        dat$p_text <- ""
        dat$p_text[1] <- p_text
        dat_tmp <- rbind(dat_tmp , dat)
    }

    y_max <- max(dat_tmp$BurdenExon) + 2
    dat_tmp$Class <- factor( dat_tmp$Class , levels = c("IGC" , "DGC") , order = T )
    #dat_tmp$From <- paste0( dat_tmp$From , "\n" , "(" , dat_tmp$nums , ")" )

    plot <- ggplot(dat_tmp, aes(x = From , y = BurdenExon, color = Class)) +
            geom_boxplot(size = 1.2 , outlier.alpha=0) + ## 去除散点，加粗线
            geom_jitter(position = position_jitterdodge(0.2) , size = 1) + 
            scale_y_sqrt() +
            scale_color_manual(values=col[c("IGC" , "DGC")]) +
            geom_text(aes(label=p_text , y = y_max , x = From),parse = TRUE,size=4 , color = "black") +
            xlab(NULL) +
            ylab(y_lab)+
            theme_bw() +
            theme(
                legend.position = 'top',
                legend.title = element_blank() ,
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 12,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold') ,
                axis.ticks.length = unit(0.2, "cm") ,
                axis.line = element_line(size = 0.5)) 
        
    out_name <- paste0(out_path , "/MutationBurden.compare.IGC_DGC.",msi,".pdf")
    ggsave(file=out_name,plot=plot,width=8/1.2,height=6/1.5)
}

###########################################################################################
## 合并总的突变负荷的比较
dat_combine$From <- "All"

for( msi in unique(unique(dat_combine$MS_Type)) ){

    dat_plot_tmp_use <- subset( dat_combine , MS_Type == msi )
    ###########################################################################################
    ## 样本数量
    sample_num <- dat_plot_tmp_use %>%
    group_by(From) %>%
    summarize( nums = length(unique(Tumor)) )

    sample_class_num <- dat_plot_tmp_use %>%
    group_by(Class) %>%
    summarize( nums_class = length(unique(Tumor)) )

    dat_plot_tmp_use <- merge( dat_plot_tmp_use , sample_num , by = "From")
    dat_plot_tmp_use <- merge( dat_plot_tmp_use , sample_class_num , by = "Class")


    dat_tmp <- c()

    for( from in unique(dat_plot_tmp_use$From) ){

        dat <- subset( dat_plot_tmp_use , From == from )

        a <- dat[dat$Class==unique(dat$Class)[1],"BurdenExon"]
        b <- dat[dat$Class==unique(dat$Class)[2],"BurdenExon"]
        
        p <- wilcox.test( a , b )$p.value

        if( p < 0.01 ){
            p_text <- trans(p)
        }else{
            p_text <- paste0( "P == " , round(as.numeric(p) , 3) ) 
        }
        dat$p_text <- ""
        dat$p_text[1] <- p_text
        dat_tmp <- rbind(dat_tmp , dat)
    }

    dat_tmp$Class_num <- paste0( dat_tmp$Class , "\n" , "(" , dat_tmp$nums_class , ")" )
    y_max <- max(dat_tmp$BurdenExon) + 12
    dat_tmp$Class_num <- factor( dat_tmp$Class_num , levels = unique(dat_tmp$Class_num)[2:1] , order = T )
    dat_tmp$Class <- factor( dat_tmp$Class , levels = unique(dat_tmp$Class)[2:1] , order = T )

    plot <- ggplot(dat_tmp, aes(x = Class_num , y = BurdenExon, color = Class , fill = Class)) +
            #geom_boxplot(size = 1.2 , outlier.alpha=0) + ## 去除散点，加粗线
            #geom_jitter(position = position_jitterdodge(1) , size = 1) + 
            geom_violin(trim=FALSE) +
            geom_boxplot(width=0.2,position=position_dodge(0.9),fill="white",color="black")+ #绘制箱线图
            #scale_y_sqrt() +
            scale_y_continuous(
                trans = sqrt_trans(),
                breaks = c(1,2,4,8,16,32)
                ) +
            scale_color_manual(values=col[c("IGC" , "DGC")]) +
            scale_fill_manual(values=col[c("IGC" , "DGC")]) +
            geom_text(aes(label=p_text , y = y_max , x = 1.5),parse = TRUE,size=5 , color = "black" , face='bold') +
            xlab(NULL) +
            ylab(y_lab)+
            theme_bw() +
            theme(
                legend.position = 'none',
                legend.title = element_blank() ,
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 12,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                axis.text.x = element_text(size = 12,color="black",face='bold') ,
                axis.ticks.length = unit(0.2, "cm") ,
                axis.line = element_line(size = 0.5)) 
        
    out_name <- paste0(out_path , "/MutationBurden.combine.IGC_DGC.",msi,".pdf")
    ggsave(file=out_name,plot=plot,width=2.3,height=3.6)

    outcompare <- dat_combine %>%
    group_by( Class , From ) %>%
    summarize( BurdenExon = median(BurdenExon) )
    out_name <- paste0(out_path , "/MutationBurden.combine.IGC_DGC.",msi,".tsv")
    write.table( outcompare , out_name , row.names = F , quote = F , sep = "\t" )
}